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Abstract 


In  this  research,  a  finite  element  method,  originally  developed  for  analyzing  the 
static  behavior  of  composite  flat  plates,  was  enhanced  so  it  could  also  be  used  with 
sandwich  plates.  The  governing  theory  considers  geometric  nonlinearity  and  transverse 
shear  effects.  Furthermore,  a  new  external  postprocessor  was  written  in  order  to  check 
plate  models  for  initial  failure  using  the  maximum  stress  criteria.  It  also  includes  a 
procedure  for  evaluating  transverse  normal  stresses  by  enforcing  equilibrium  through  the 
thickness. 

The  programming  modifications  to  allow  modeling  of  sandwich  plates  were 
verified  by  comparing  finite  element  solutions  to  those  from  closed-form  linear  theories 
for  sandwich  plates.  Results  showed  good  correlation  between  the  numerical  and 
theoretical  solutions.  In  addition,  displacement  results  (using  the  same  program)  from 
previous  research  for  particular  composite  plates  were  compared  to  sandwich  plates  of 
similar  composition  and  equal  size.  The  sandwiches  were  more  flexible  in  absolute  terms, 
but  displayed  higher  stiffness-to-weight  ratios  than  the  composites.  Finally,  low-velocity 
impact  tests  were  modeled  quasi-statically  with  the  finite  element  code,  and  the  new 
postprocessor  was  employed  to  predict  incipient  plate  damage.  Locations  and  modes  of 
failure  were  correctly  determined,  but  the  predicted  load  levels  for  initial  failure  were 
inconsistent  with  experimental  results  from  other  research. 


FTNTTF.  F.T.F.MF.'MT  methods  for  nonlinear  static  analysis  of  sandwich  plates 


I.  Introduction 


The  need  for  strong,  lightweight  materials  in  aerospace  structural  components  has 
led  to  a  renewed  interest  in  sandwich  plates  and  shells  using  modem  composite  materials. 
These  hybrid  constructions  consist  of  two  dense  outer  faces  that  are  bonded  to  a 
lightweight  core.  The  core  usually  has  little  in-plane  and  flexural  stiffness,  compared  to 
the  faces,  but  it  can  have  significant  transverse  strength  and  acts  as  a  spacer  to  enhance 
the  bending  resistance  of  the  faces.  The  result  is  a  thicker  plate  or  shell  with  a  higher 
stiffness-to-weight  ratio  than  the  facesheets  alone. 

Predicting  the  static  response  of  laminated  composite  plates  is  complicated  due  to 
effects  such  as:  property  variation  through-the-thickness,  geometric  and  physical 
nonlinearity,  transverse  shear  and  multiple  failure  modes.  The  additional  complexity  of 
sandwich  constructions  further  complicates  the  analysis.  Closed-form  methods  are 
limited  to  linear  solutions  (with  many  simplifying  assumptions)  for  specific  geometries 
and  boundary  conditions.  Furthermore,  experimental  testing  can  give  good  results  for  a 
particular  plate,  but  it  can  be  impractical,  in  terms  of  time  and  money,  for  analyzing  the 
effects  of  a  wide  range  of  variables.  On  the  other  hand,  numerical  techniques  like  finite 
elements  (FE)  can  be  applied  to  plates  of  different  shapes,  sizes,  compositions,  loadings 
and  supports  with  greater  flexibility.  The  accuracy  and  practicality  of  FE  methods  are 
dependent  on  the  governing  theories,  model  complexity  and  a  given  computer  s  speed 
and  precision.  Therefore,  the  application  of  advanced  theories  and  finely-detailed  models 
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is  conditioned  by  the  availability  of  affordable  technology  that  can  handle  them  in  a 
timely  manner. 

This  research  employs  a  finite  element  method  developed  for  laminated  plates  and 
cylindrical  shells  and  enhances  it  for  use  with  sandwich  plates.  The  governing  theory 
considers  geometric  nonlinearity  (through  the  von  Karman  strain-displacement  relations) 
and  transverse  shear  effects.  In  addition,  plate  materials  are  assumed  linearly  elastic  to 
exclude  such  effects  as  plasticity  and  failure  within  the  solution.  A  full  three-dimensional 
plate  model  is  reduced  to  a  two-dimensional  analysis  by  describing  all  displacement 
variations  in  the  thickness  direction  relative  to  those  at  the  mid-plate  surface.  This 
assumption  ignores  transverse  normal  stresses,  but  it  greatly  reduces  the  solution’s 
complexity. 

The  main  objective  of  this  research  was  to  test  the  effectiveness  of  the  given  finite 
element  method  in  analyzing  sandwich  plates.  Three  cases  were  considered.  First,  linear 
and  nonlinear  displacement  results  were  obtained  for  the  same  sandwich  plate  models 
used  by  Pagano  [18]  and  Whitney  [21]  in  developing  closed-form  solutions.  This 
provided  comparisons  to  established  theories  in  order  to  verify  the  FE  algorithms. 

Second,  displacement  results  (using  the  same  FE  code)  from  previous  research  by  Owens 
[16,17]  for  laminated  plates  were  compared  to  sandwiches  of  equivalent  overall 
geometry-  in  which  the  core  was  half  the  total  thickness  and  the  face  plies  were 
constructed  from  the  same  material  as  the  laminated  plates.  Finally,  an  attempt  was  made 
to  predict  initial  failure  using  postprocessed  stress  calculations  and  maximum  stress 
failure  criteria.  Sandwich  plates  used  in  experimental  work  by  Harrington  [7]  were 
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modeled,  so  that  the  FE  results  could  be  compared  to  actual,  incipient  plate  damage  from 
low  velocity  impacts.  In  this  research,  a  quasi-static  approach  was  used  to  simulate  the 
dynamic  loading  employed  in  the  experiments. 

Previous  Work 

The  textbook  by  Palazotto  and  Dennis  [19]  gives  a  detailed  history  of  the  use  of 
finite  elements  for  analyzing  flat  plates  and  cylindrical  shells.  It  also  describes  the  past 
research  and  theories  which  led  to  the  development  of  the  FE  code  (called  SHELL)  used 
in  this  thesis.  Dennis  [4]  wrote  the  original  version  of  SHELL  for  the  study  of  large 
displacements  and  rotations  of  shells.  Owens  [16]  used  it  to  analyze  composite  plates  and 
made  comparisons  between  linear,  geometrically  nonlinear,  classical  and  nonclassical 
solutions.  He  showed,  for  a  given  loading,  that  membrane  stiffness  due  to  nonlinear  in¬ 
plane  strain  terms  becomes  significant  in  thin  plates.  In  addition,  the  transverse  shear 
flexibility  present  in  nonclassical  theories  is  important  for  thick  plates.  Linear  and 
nonlinear  solutions  become  alike  for  thicker  plates  as  do  classical  and  nonclassical 
solutions  for  thinner  plates. 

Early  work  in  the  study  of  sandwich  plates  was  conducted  by  Pagano  [1 8],  He 
developed  a  linear,  three-dimensional  elasticity  solution  for  rectangular  plates  with 
simply-supported  edges.  Pagano’s  results  for  both  composite  and  sandwich  plates  further 
emphasized  the  limitations  of  classical  laminated  plate  theory  (CLPT)  for  thick  plates  due 
to  its  neglect  of  transverse  shear.  Later  work  by  Whitney  [21]  yielded  an  alternative 
closed-form  solution  that  resembles  CLPT  with  additional  contributions  for  transverse 
shear.  Displacement  results  from  both  methods  showed  good  agreement  over  a  wide 
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range  of  plate  width-to-thickness  ratios;  although,  for  thin  plates,  Pagano’s  solution 
converged  closer  to  CLPT  than  Whitney’s  did. 

Both  Pagano  and  Whitney  calculated  in-plane  stresses  from  strain  values  and  the 
constitutive  relations,  but  they  derived  transverse  shear  stresses  from  the  equilibrium 
equations  in  elasticity  theory  (see  the  textbook  by  Sadda  [20]).  This  method  provides 
better  accuracy  than  computing  all  stresses  through  the  constitutive  relations.  Engblom 
and  Ochoa  [5,6]  also  used  this  procedure  in  developing  linear  finite  element  formulations 
for  laminated  composite  plates.  By  assuming  linear  stress  distributions  through  the 
thickness  of  each  ply,  the  integral-differential  equations  of  equilibrium  can  be  converted 
into  matrix  operations  for  easy  implementation. 

On  the  other  hand,  SHELL  does  not  satisfy  localized  equilibrium  (although  such 
errors  tend  to  cancel  out  on  a  global  scale)  since  it  relates  all  stresses  to  strains  through 
the  constitutive  relations.  It  does  this  because  its  consideration  of  nonlinear  strains  and 
higher-order  stress  distributions  prevents  it  from  obtaining  transverse  stresses  within  the 
solution  through  a  linear  system  of  equations  that  enforce  equilibrium.  Furthermore, 
compatibility  is  satisfied  at  the  nodes  but  not  through  the  thickness,  in  general.  SHELL 
could  be  modified  to  employ  an  iterative  process  in  which  the  postprocessor  calculates 
transverse  stresses  by  enforcing  equilibrium  and  then  uses  them  to  alter  the  strain  and 
displacement  solutions  until  both  equilibrium  and  compatibility  are  satisfied.  However, 
this  would  have  to  be  nested  within  the  existing  iteration  scheme  for  load  or  displacement 
incrementing  in  SHELL’S  nonlinear  solution  control,  thus  greatly  reducing  the  code’s 
speed.  Furthermore,  such  a  modification  is  beyond  the  scope  of  this  thesis. 
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Other  research  related  to  sandwich  plates  has  been  oriented  towards  studying  and 
predicting  failure  due  to  damage  brought  about  by  low  velocity  impacts.  McQuillan  et  al. 
[14]  were  among  the  first  to  observe  that  static  and  dynamic  loading  had  similar  failure 
mechanisms.  Experimental  work  by  Kelkar  et  al.  [10]  employed  this  concept  to  simulate 
low-velocity  impact  damage  using  equivalent  quasi-static  loads.  Low-velocity  impacts 
are  representative  of  physical  damage  caused  by  such  actions  as  a  dropping  a  blunt  object 
onto  the  plate  from  a  short  height.  In  addition,  quasi-static  means  that  loads  are  applied 
slowly  enough  to  ignore  inertial  effects,  and  the  plate  is  assumed  to  remain  in  static 
equilibrium  as  damage  progression  alters  its  equilibrium  state.  Finite  element  models 
using  ANSYS  (a  commercial  software  package)  were  developed  by  Dandy  et  al.  [3]  for 
comparison  with  Kelkar’s  experiments.  Both  had  agreeable  results  for  thick  plates  but  not 
for  thin  plates.  The  ANSYS  solver  can  only  consider  linear  strains,  and  the  large 
displacements  present  in  the  thin  plates  invalidate  this  simplification. 

Nemes  and  Simmonds  [15]  analyzed  the  impact  response  of  sandwich  plates  with 
a  foam  core.  The  experimental  plates  were  square,  but  the  finite  element  models  were 
circular  disks  of  equivalent  size.  This  allowed  a  complicated  three-dimensional  model  to 
be  reduced  to  an  axisymmetric  radial  plane-section  of  the  disk.  Therefore,  each 
sandwich’s  cross-section  through  the  thickness  was  modeled  as  a  continuum  instead  of 
distinct  layers.  The  FE  results  overpredicted  the  experimental  deflection  but  produced 
agreeable  transverse  shear  stresses.  In  addition,  the  low  stiffness  of  the  foam  cores 
allowed  significant  relative  displacement  between  the  faces.  This  effect  could  be  reduced 
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by  introducing  a  honeycomb  core  with  greater  transverse  normal  stiffness,  but  the  voids 
present  in  honeycomb  cells  may  invalidate  treating  it  as  a  continuum. 

Sandwich  plate  failure  modes  due  to  impacts  and  static  indentations  were 
investigated  by  Lagace  and  Williamson  [12].  As  predicted,  both  types  of  loadings  had 
similar  damage  responses.  They  tested  sandwiches  with  graphite-epoxy  laminate  faces 
and  Nomex  honeycomb  cores,  which  are  similar  to  those  used  in  this  thesis.  Initial 
damage  occurred  under  the  applied  load  from  core  crushing  or  buckling  near  the  top  face. 
As  a  consequence,  this  face  experienced  local  fiber,  matrix  and  delamination  failure 
because  the  core  could  no  longer  support  it.  Thicker  faces  reduced  the  extent  of  damage 
by  stiffening  the  indentation  responses.  On  the  other  hand,  thicker  cores  increased  the 
chances  of  buckling  more  than  they  enhanced  the  faces’  bending  stiffness. 

In  summary,  the  majority  of  research  dealing  with  sandwich  plates  has  been 
limited  to  experimentation  and  linearized  solutions.  This  thesis  goes  one  step  further  by 
employing  a  finite  element  method  that  includes  geometric  nonlinearities.  The  goal  is  to 
validate  it  as  an  accurate  and  practical  tool  for  static  displacement  and  initial  failure 
analysis. 
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II.  Theoretical  Considerations 


Sandwich  Plate 

Figure  2.1  contains  the  geometry  and  coordinate  systems  used  for  modeling 
sandwich  plates.  Both  X-Y-Z  (XrX2-X3)  and  L-T-Z  represent  orthogonal  systems.  The 
longitudinal  and  lateral  directions  correspond  to  the  principal  material  directions  of  an 
orthotropic  ply.  A  ply’s  orientation  angle  0  is  the  angle  from  X  to  L  (or  from  Y  to  T).  All 
plates  used  in  this  research  were  symmetrical  about  their  midsurfaces  (z=0). 


Figure  2.1:  Sandwich  Plate  Geometry  and  Coordinate  Systems 
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Governing  Equations 

A  plate  is  assumed  to  be  in  a  state  of  plane  stress.  As  a  result,  all  transverse 
normal  stresses  ctzz  are  zero,  and  plate  behavior  can  be  described  by  displacements  and 
rotations  at  and  relative  to  the  midsurface.  Transverse  normal  strains  8ZZ  are  nonzero  in 
general,  but  they  are  consequences  (due  to  Poisson  effects)  of  the  other  strains  and  do  not 
affect  the  stress  state.  Transverse  shear  strains  sxz  and  syz  are  assumed  to  have  parabolic 
distributions  in  the  Z-direction.  This  can  be  done  since  the  plane  stress  assumption 
decouples  the  in-plane  and  transverse  shear  constitutive  relations.  It  also  satisfies  the 
boundary  conditions  of  zero  transverse  shear  on  the  top  and  bottom  plate  surfaces  (none 
of  the  prescribed  loading  in  this  research  imposed  surface  shears). 


Each  point  within  the  plate’s  midsurface  has  seven  degrees-of-freedom  as  shown 
in  Figure  2.2.  Displacements  u,  v  and  w  are  translations  in  the  X,  Y  and  Z  directions.  The 
terms  w5l  and  w,2  are  physical  slopes  of  the  midsurface  in  the  X-Z  and  Y-Z  planes,  while 
vg!  and  V) ;2  are  rotations  due  to  bending  alone  in  those  respective  planes.  Transverse 
shearing  in  a  single  plane  is  described  by  the  algebraic  sum  of  the  two  rotations. 
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Translational  displacements  away  from  the  midsurface  are  evaluated  through  the 
following  plate  kinematics: 


u{(x,y,z)-u  +  z\\i  {  +z3£(y{/]  +wh  ) 
u2(x,y,z)  =  v  +  z\\f  2  +z3k(\\i  2  +w,2 ) 
u3(x,y,z)  =  w 

k--M  (3h2 ) 


Furthermore,  nonlinear  strain  and  displacement  are  related  through  the  von  Karman  plate 
equations  [19]  (linear  plate  solutions  disregard  all  nonlinear  terms): 

eXx  =um  +  \w’i2 

£yy  =U2’ 2+2  W>22 

8  Xy  —  Wj  ,2  +W2  ?i  +W5]  XV, 2  (2-2) 

£yz=(l  +  3z2k)(w,2+\v  2) 

Exz  =(1  +  3z1k,)(w,x+\y  1) 

All  ply  materials  are  assumed  linearly  elastic  and  at  least  orthotropic.  The 
assumption  of  plane  stress  allows  the  need  for  only  six  elastic  constants:  ELL  ,  EXT  ,  GLT  , 
Gtz  5  Glz  and  vLX.  The  constitutive  relations  for  stress  and  strain  are: 
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where  K  is  the  ply  number  and  Q?  are  the  components  of  that  ply  material’s  elastic 
stiffness  matrix--  reduced  for  plane  stress  and  transformed  into  the  X  and  Y  directions. 
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The  numerical  subscripts  represent  a  simplified  indexing  of  the  stress  and  strain  tensor 


components: 
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Finite  Element  Solution 

This  research  employed  rectangular  plate  elements  with  four  nodes  and  28 
degrees-of-freedom  (seven  per  node  as  in  Figure  2.2).  The  geometry  of  an  individual 
element  and  the  representation  of  its  global,  local  and  natural  coordinates  are  shown  in 
Figure  2.3.  Displacements  within  the  given  element  are  interpolated  from  the  nodal 
displacements  through  appropriate  shape  functions.  The  displacement  field  for  w 


Local 


X 

/N 


Natural 


X 

/N 


Ax 


Global 


N'i 


(Ax/2 ,  Ay/2) 


-»Y 


Node  2 

Node  3' 

Node  1 

Node  4 

.  4 

/ 

\ 

\ 

Ay  / 

Figure  2.3:  Four-Node  Plate  Element  Geometry  and  Coordinate  Systems 
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requires  C1  continuity  (as  defined  in  the  textbook  by  Cook  et  al.  [2]),  therefore  Hermitian 
shape  functions  are  used  for  nodal  displacements  w  ,  w5l  and  w,2: 


w(x,y)=[H,  H2  H 3  HA 


<h 
<h 
<h 
1*4  J 


^K  = 


(l+^K^)(l+riKTi)(2+^K^+TiKTi-^-n2) 

iAx^K(l+^K02(^^-l)(l+^l^) 

l  Ayr)  K(l+£,  k£,  )(rj  Kty -l)(l+r|  Kr)  )2 
qK  =  {w  w„  w,2}* 


(2.5) 


where  K=1  through  4  represent  the  local  node  numbers  for  an  element  found  at  global 
position  (x,y).  The  other  displacement  fields  only  need  C°  continuity  and  employ 
Lagrangian  shape  functions: 
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The  complex  formulation  of  SHELL’S  finite  element  solution  is  fully  described  in 
the  textbook  by  Palazotto  and  Dennis  [19].  In  a  nonlinear  model,  the  code  allows  either 
load  or  displacement  incrementing  for  solution  control.  For  each  increment,  it  uses  a 
Newton-Raphson  iteration  scheme  to  converge  to  a  solution  which  minimizes  potential 
energy.  All  plates  modeled  in  this  research  used  some  form  of  distributed  or  multi-nodal 
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loading.  Therefore,  load  control  was  employed  in  each  nonlinear  case,  since 
displacement  control  would  require  an  assumed  shape  relation  between  the  prescribed 
nonzero  degrees-of-freedom.  In  this  research,  no  iteration  convergence  problems  arose 
using  load-control. 

Every  case  study  in  this  research  considered  square  plates  with  simply  supported 
edges  (u  and  v  translations  were  free).  Since  all  ply  orientations  were  either  0  or  90 
degrees,  it  was  only  necessary  to  generate  FE  meshes  for  a  single  quadrant  of  each  plate 
by  prescribing  bi-axial  symmetry.  Figure  2.4  shows  the  displacement  boundary 
conditions  that  were  applied  to  each  square  quarter-plate. 
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Figure  2.4:  Boundary  Conditions  of  Square  Quarter  Plate 
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Code  Enhancements 


SHELL  required  some  modifications  before  being  used  with  sandwich  plates.  The 
most  critical  change  was  to  allow  multiple  sets  of  elastic  properties.  Without  this,  the 
faces  and  core  could  not  be  represented  as  different  materials.  Furthermore,  the  typically 
large  variation  between  face  and  core  thickness  commanded  the  need  to  define  a  separate 
thickness  for  each  ply.  Otherwise,  it  would  be  necessary  to  divide  all  plies  into  a  common 
uniform  thickness,  which  could  cause  severe  redundancies.  Finally,  SHELL  now  gives 
the  user  the  opportunity  to  generate  a  secondary  output  file  for  use  with  a  newly  written 
(as  part  of  this  research)  and  separately  executed  postprocessor  program,  called 
FAILURE.  Appendix  A  provides  a  more  detailed  explanation  of  these  and  other 
enhancements  to  SHELL  and  includes  the  new  structure  of  its  input  deck. 

FAILURE  was  written  for  the  purpose  of  predicting  the  initial  failure  regions  and 
modes  of  rectangular  plates  (with  certain  modeling  restrictions)  using  the  maximum 
stress  criteria.  SHELL’S  own  postprocessor  could  have  been  altered  for  this  task.  In  fact, 
FAILURE  utilizes  some  of  the  same  subroutines  (with  minor  modifications).  However,  a 
separate  program  has  several  advantages.  First,  the  same  plate  model  can  be  rechecked 
for  failure  using  different  parameters  and  criteria  values  without  having  to  re-execute 
SHELL  to  obtain  the  same  displacement  solution  before  postprocessing  (a  valuable 
feature  for  code  debugging  and  validation  of  the  methodology,  since  it  has  not  been  tried 
before  with  this  FE  theory).  Second,  the  code’s  structure  does  not  need  to  conform  to  that 
of  SHELL.  Therefore  it  gives  the  user  more  flexibility  in  adding  or  changing  program 
features.  On  the  other  hand,  the  secondary  output  file  generated  by  SHELL  (which 
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contains  all  preprocessor  and  solution  data  needed  by  FAILURE)  is  usually  as  large  as  its 
regular  output  file.  Hence,  using  it  with  a  set  of  complex  FE  models  requires  a  computer 
with  plentiful  storage  space.  Appendix  B  includes  the  entire  FAILURE  code,  the 
structure  of  its  additional  user-defined  input  deck  and  some  of  its  other  features. 

Since  postprocessor  results  are  based  on  the  assumption  of  perfect  linearly  elastic 
materials,  FAILURE  is  invalid  for  predicting  failure  beyond  the  point  of  initiation  (even 
with  a  more  sophisticated  failure  theory  than  maximum  stress).  In  addition,  SHELL’S 
basic  design  does  not  allow  elastic  property  variation  from  element-to-element  (a 
necessary  feature  for  localizing  the  effects  of  failure  within  the  solution).  Hence,  the 
implementation  of  progressive  failure  into  the  FE  solution  would  require  massive 
amounts  of  code  alteration. 


Initial  Failure  Criteria 


FAILURE  is  designed  to  report  failure  when  averaged  stresses  within  a  given 
element  exceed  user-defined  maximum  magnitudes.  For  a  laminated  plate,  element 
stresses  are  calculated  at  12  discrete  points  per  ply-  the  four  outermost  Gauss  points  (see 


Figure  2.5)  at  the  upper,  middle  and  lower  surfaces  of  each  ply  (the  FE  solution  uses  5x5 
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Note:  Outer  Gauss  points  are 
located  near  comer  nodes  but  are 
not  shown  to  scale.  See  [2]  for 
information  on  numerical 
positioning. 


Figure  2.5:  Single-Element  Outer  Gauss  Points  for  Stress  Calculation 
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Gauss  quadrature  for  numerical  integration  [2,19]).  Therefore,  at  a  given  Gauss  point, 
each  ply’s  stress  distribution  through  the  thickness  is  characterized  by  three  values  per 
component  at  known  Z-coordinates.  This  provides  sufficient  data  to  obtain  an  average 
stress  state  by  calculating  parabolas  that  fit  each  component’s  discrete  quantities  and  then 
finding  a  mean  value  for  each  continuous  function  (the  area  under  the  curve  divided  by 
the  ply  thickness).  This  is  a  better  method  because  it  is  equivalent  (for  the  assumed  shape) 
to  the  arithmetic  average  of  an  infinite  number  stress  values  (per  component)  instead  of 
just  three.  Figure  2.6  graphically  demonstrates  how  ply  stresses  are  averaged. 


°avg 


Figure  2.6:  Ply  Stress  Averaging  Through  the  Thickness 


Average  stress  components  are  given  in  X-Y-Z  coordinates,  so  they  must  be 
transformed  to  match  the  directions  of  the  material  failure  criteria.  For  an  orthotropic  ply, 
the  stresses  are  transformed  into  the  L-T-Z  system  by  the  following  matrix  operation: 
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where  c=cos  0  and  s=sin  0.  On  the  other  hand,  stresses  for  an  isotropic  ply  are  converted 
to  principal  stresses  by  evaluating  the  eigenvalues  of  the  stress  tensor-  as  shown  in  the 
textbook  by  Sadda  [20].  The  resulting  characteristic  polynomial  is  cubic  and  its  roots  are 
solved  through  a  closed-form  technique  found  in  the  mathematics  handbook  by  Korn  and 
Korn  [11].  Seven  possible  modes  of  failure  were  given  to  orthotropic  materials  and  three 
were  given  to  isotropic  materials.  These  modes  and  their  maximum  stress  criteria  are 
listed  in  Table  2.1. 


Table  2.1 :  Material  Failure  Modes  and  Criteria 


Material 

Mode 

Criteria 

Orthotropic  (L-T-Z  stresses) 

Longitudinal  Tension 
Longitudinal  Compression 

Lateral  Tension 

Lateral  Compression 
Long.-Lat.  Shear 

Long.-Z  Shear 

Lat.-Z  Shear 

*TlL  —  ^LL  max 

aLL  -  °LL  min 

aTT  >  axx  max 

Gxx  <  G  rx  min 

I^LtI  -°LTmax 
laLzl  -  aLZ  max 
|<7Xzl  -  aTZ  max 

Isotropic  (principal  stresses) 

Tension 

Compression 

Shear 

°p  max  —  ^max  uniaxial 

°p  min  —  '^max  uniaxial 

max  "  min  —  ®max  uniaxial 

FAILURE  also  includes  a  procedure  for  checking  delamination  at  the  ply 
interfaces  due  to  transverse  shearing.  Since  the  constitutive  relations  (Equation  2.3) 
cause  stress  discontinuities,  in  general,  across  the  interface  between  unlike  plies,  the 
values  of  axz  calculated  on  each  side  of  the  interface  are  usually  different.  This  also 
holds  true  for  crw.  The  arithmetic  mean  stress  at  the  interface  for  each  component  is  used 
for  checking  delamination.  FAILURE  reports  shear  delamination  for  a  given  element  and 
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interface  if  (at  any  Gauss  point)  the  magnitude  of  either  mean  value  exceeds  defined 
maximums. 

Estimation  of  Transverse  Normal  Stresses 

Although  SHELL’S  theory  assumes  all  transverse  normal  stresses  are  zero, 
FAILURE  includes  a  routine  which  estimates  azz  by  enforcing  equilibrium.  Neglecting 
body  forces,  the  equation  of  equilibrium  in  the  Z-direction  is: 

°xz>x+CTyz>y+CTzz>z  =  0  (2-8) 

The  shear  stress  gradients  in  Equation  2.8  are  related  to  displacement  gradients  of  w , 

\j/j  and  \|/2  by  taking  partial  derivatives  of  Equations  2.2  and  2.3: 
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(2.10) 


The  displacement  gradients  in  Equation  2.10  are  related  to  nodal  displacements  through 
derivatives  of  the  shape  functions  in  Equations  2.5  and  2.6.  Since  these  gradients  are  also 
used  in  other  strain  terms,  SHELL’S  stress  calculation  subroutine  (which  is  modified  for 
use  with  FAILURE)  already  contains  code  for  determining  their  values.  Furthermore,  the 
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assumed  parabolic  distribution  of  transverse  shear  strains  forces  each  ply’s  transverse 
shear  stresses,  and  their  in-plane  gradients,  to  be  parabolic  functions  of  Z: 

o  yZjy  (x,  y,  z )  =kA(x,  y)z 2  +K5(x,  y)z+KC(x ,  y) 
xz,x(^>  T*  z)=kD(x,  y)z2  +kE(x,  y)z+KF(x,  y) 

where  z  is  located  within  a  given  ply  K.  The  parabola  coefficients  for  a  particular  ply  and 
Gauss  point  are  obtained  by  curve  fitting  the  discrete  values  of  ayz,y  and  oxz,x  at  the  ply’s 
upper,  middle  and  lower  ply  surfaces. 

Equation  2.8  forces  ctzz,z  =  0  at  a  plate’s  top  and  bottom  surfaces,  z  =  ±h/2, 
because  ayz  (x,  y,  +h/2)  =  axz  (x,  y,  ±h/2)  =  0  from  Equations  2.2  and  2.3.  In  other  words, 
since  the  transverse  shear  stresses  are  zero  everywhere  on  those  surfaces,  their  in-plane 
gradients  must  also  be  zero  on  those  surfaces.  In  addition,  the  finite  element  solution 
assumes  that  all  prescribed  transverse  loading  occurs  on  the  top  surface.  Therefore,  the 
bottom  of  the  plate  is  free  of  transverse  normal  stresses:  azz  (x,  y,  h/2)  =  0.  By  combining 
Equations  2.8  and  2.1 1  and  integrating  in  the  Z-direction,  the  calculation  of  azz  becomes: 


=  J[(K^+Kfl)z!  +eB+KE)z+(KC+*F)]dz 

K=I  *K-1 

=  X  ti (kA+kD)z3  +\(kB+kE)z2  +(kC+kF)z] 


(2.12) 


K=1 


ZK 

ZK-1 


where  K=1  denotes  the  top  ply  and  K=n  is  the  ply  located  at  z=  -zn.  The  bounds  zKA  and 
zK  are  the  Z-coordinates  of  ply  K’s  upper  and  lower  surfaces,  except  when  K=n.  In  that 
case  zK=  -zn.  Note  that  the  use  of  -zn  arises  because  the  Z-axis  is  positive  downward 
through  the  plate.  Integrating  in  the  positive  Z-direction  imposes  an  initial  value  of  gzz=0 
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at  the  top  of  the  plate  (z=  -h/2).  Equation  2.12  includes  sign  changes  that  alter  the 
constant  of  integration  (and  thus  translates  the  function  gzz(z)  for  a  fixed  x  and  y)  in  order 
to  enforce  the  boundary  condition  of  ctzz=0  at  the  bottom  of  the  plate  (z=  h/2).  Code 
testing  verified  that  Equation  2.12  calculates  negative  (compression)  values  of  azz(x, 
y,-h/2)  in  plate  regions  subjected  to  a  compressive  transverse  pressure  on  the  top  surface. 

Plate  symmetry,  with  respect  to  the  midsurface,  forces  all  transverse  shear  stress 
profiles  in  the  Z-direction  (and  their  in-plane  gradients)  to  be  symmetric  about  z=0. 
Hence,  integrating  Equation  2.8  generates  crzz  profiles,  for  a  fixed  x  and  y,  that  are  the 
superposition  of  an  antisymmetric  function  of  z  and  a  constant  (a  line  normal  to  z=0).  As 
a  consequence,  the  three  czz  boundary  conditions  on  the  top  and  bottom  plate  surfaces 

(czz,z(x,y,±h/2)=azz(x,y,h/2)=0)  cause  each  profile  to  resemble  a  cubic  polynomial  with  a 
maximum  magnitude  at  the  top  surface  of  ozz(x,y,-h/2).  However,  for  a  laminate,  each 
ply  is  usually  associated  with  a  different  cubic  polynomial  because  the  coefficients  in 
Equation  2.1 1  change  from  ply  to  ply.  Therefore,  the  actual  profiles  generated  from 
Equation  2.12  are  continuous  but  piecewise  smooth  at  the  ply  interfaces.  Figure  2.7 
displays  a  typical  ctzz  profile  for  a  sandwich  plate  region  subjected  to  a  transverse 
pressure  on  the  top  surface. 

Testing  of  the  FAILURE  code  revealed  a  numerical  problem  associated  with 
directly  calculating  cryz,y  and  crxz,x  from  nodal  displacements.  The  use  of  Lagrangian 
shape  functions—  with  only  C°  continuity—  for  u  ,  v  ,  i|/ j  and  v|/2  (Equation  2.6)  cause  in¬ 
plane  discontinuities  in  all  stress  fields.  The  effects  are  usually  less  significant  for  oxx , 
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oyy  and  axy  because  their  corresponding  strains  have  more  w  ,  wsl  and  w,2  terms 
(Equations  2.1  and  2.2)  which  tend  to  minimize  the  variations  across  adjacent  elements. 


Figure  2.7:  Shape  of  a ^  Distribution  for  Sandwich  Plate  under  Transverse  Pressure 


However,  transverse  shear  strains  can  fluctuate  severely  in  the  direction  normal  to  the 
transverse  plane-  the  X-direction  for  syz  and  the  Y-direction  for  sxz.  Therefore,  the 
fluctuations  amplify  when  calculating  cyz,y  and  gxz,x  from  Equation  2.9  because  they 
require  transverse  strain  gradients  in  both  X  and  Y-directions.  This  resulted  in  highly 
varied  and  inaccurate  estimations  of  ctzz.  Fortunately,  a  simple  averaging  of  the  ayz,y  and 
oxz,x  values  obtained  at  an  element’s  outer  Gauss  points  greatly  reduced  these 
fluctuations.  Using  the  Gauss  point  labeling  as  shown  in  Figure  2.5,  the  eight  stress 


2-14 


gradient  values  for  a  given  element  and  ply  surface  (upper,  middle  or  lower)  are 
combined  into  four  average  values. 
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(2.13) 


When  the  calculated  stress  gradients  were  replaced  by  these  averages,  sample  plates  with 
a  uniform  transverse  pressure  obtained  peak  crzz  values  (from  Equation  2.12)  that 
typically  ranged  between  65  and  85  percent  of  the  applied  pressure.  These  estimates  were 
better  than  expected,  since  ctzz  is  related  to  the  nodal  displacement  through  third-hand 
calculations  (strain-displacement  equations,  constitutive  relations  and  equilibrium  in  the 
Z-direction). 

One  problem  with  Equation  2.12  is  that  it  prevents  satisfaction  of  a  fourth 
boundary  condition,  ctzz  (x,  y,  -h/2)  =  0,  when  a  region  of  the  top  surface  is  free  of 
transverse  loading.  It  will  generally  calculate  a  nonzero  value  because  a  cubic  polynomial 
is  fully  constrained  by  four  boundary  conditions,  and  allowing  azz  (x,  y,  ±h/2)  =  ozz,z  (x, 
y,  ±h/2)  =  0  causes  zero  transverse  normal  stress  throughout  the  thickness.  The  only  way 
this  can  occur  is  if  ayz  y  (x,  y,  z)  =  -axz  x  (x,  y,  z)  for  a  particular  x  and  y. 

However,  the  previously  mentioned  research  by  Engblom  and  Ochoa  [5,6]  may 
provide  a  means  of  masking  the  problem.  They  obtained  axz  and  ayz  through  the 
equations  of  equilibrium  in  the  X  and  Y-direetions  and  also  utilized  Equation  2.8  to 
calculate  azz.  The  resulting  thickness  profile  of  ozz  for  a  stress-free  top  surface 
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resembled  a  sine  wave.  Hence,  one  possibility  is  to  map  the  curve  obtained  from 
Equation  2.12  to  a  sine  wave  with  an  equivalent  area  under  the  curve,  as  shown  in  Figure 
2.8.  The  area  calculated  from  Figure  2.7  (and  other  curves  based  on  Equation  2.12)  is 
(h/2)«G2Z(x,  y,  -h/2)  because  its  antisymmetric  shape  has  the  same  area  as  a  triangle 
formed  by  the  Z-axis,  z=  -h/2  and  the  line  connecting  cr2Z(x,  y,  h/2)  and  gzz(x,  y,  -h/2). 
Also  note  that  Figure  2.8  fails  to  satisfy  the  boundary  conditions  of  gzz,z  =  0  at  z=±h/2.  A 
more  complicated  mapping  function-  one  that  permits  additional  constraints-  could 
solve  this  problem,  but  for  simplicity  the  regions  of  small  stress  gradients  near  the  top 
and  bottom  plate  surfaces  are  assumed  negligible  with  respect  to  the  entire  thickness. 
FAILURE  does  not  presently  include  any  mapping  technique  since  the  structure  of 
SHELL’S  secondary  output  file  does  not  contain  information  on  the  location  of  nodal 
loads  used  to  generate  a  displacement  solution.  Thus,  FAILURE  cannot  tell  where 
mapping  should  be  used.  The  task  of  interpreting  the  meaning  of  the  azz  calculations  is 
left  to  the  user. 
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III.  Pagano /  Whitney  Sandwich  Plate  Models 


Closed-Form  Solutions 

Both  Pagano  [18]  and  Whitney  [21]  obtained  linear  solutions  for  simply- 
supported  square  sandwich  plates  that  are  subjected  to  sinusoidal  pressures.  Each  face 
sheet  was  a  single  [0°]  ply  of  an  unidentified  composite  material,  and  the  core  was  some 
type  of  transversely  isotropic  material.  The  thickness  of  a  single  face  was  one-tenth  that 
of  the  core  (hf  =  hc  /10  and  h=  2  hf  +  hc).  The  face  and  core  materials  had  the  following 
relevant  elastic  properties  (converted  to  SI  units): 

Face:  ELL=172.3  GPa  ETT=6.985  GPa 

GLt=Glz=3.447  GPa  GTZ=1.379GPa  (3.1) 

vLT=0.25 

Core:  EXx=EYy=275.8  MPa 

Gxy=1  10.3  MPa  Gxz=Gyz=413.7  MPa  (3.2) 

vxy==0-25 

Pagano  [18]  develops  his  elasticity  solution  for  a  generally  loaded  plate,  but  the 
results  for  these  particular  plates  are  given  numerically  at  selected  locations.  However,  he 
also  includes  a  CLPT  solution  algorithm  that  can  be  conveniently  applied  to  these  cases. 
For  the  square  quarter-plate  notation  shown  in  Figure  3.1,  the  pressure  distribution  q(x,y) 
on  the  top  surface  of  the  plate  is: 

q(x,  y)=q0  cos(ttx  /  a)  cos(rcy  /  a)  (3.3) 
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Figure  3.1:  Boundary  Conditions  of  Square  Quarter  Plate 


where  q0  is  the  peak  pressure  at  the  plate’s  center  (x=y=0).  The  transverse  displacement 
field  w(x,y)  for  a  simply  supported  plate  can  be  approximated  with  a  similar  sinusoidal 
distribution: 

w( x,  y)  =  wc  cos(ttx  /  a)  cos(7ty  /  a)  (3 .4) 

where  wc  is  the  deflection  at  the  plate’s  center.  For  a  symmetric  plate,  Pagano’s  CLPT 
solution  for  wc  becomes: 

w  = - — - r  (3.5) 

c  (Du  +  2(DU  +  2D66 )  +  D22  )(n  /  a)4 

where  Dn  ,  D22  ,  D12  and  D66  are  components  of  the  bending  stiffness  matrix  and  are 
integral  functions  of  the  constitutive  relations  in  the  Z-direction  (as  defined  in  the 
textbook  by  Agarwal  and  Broutman  [1]). 
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Whitney’s  plate  solution  for  the  same  load  and  displacement  distributions  as 


Equations  3.3  and  3.4  is  the  following  expression  [21]: 


w 


C 


g0a2 

ti2D 


(AA-A) 


D  =  AXA4A6  +  2 A2A3A5  -  AXA2  -  A4A32  -  A6A2 
A  =  Dn  +  D66  +  Gxzcmeha2  /  n 2 

A  =  Dn  +  D66 

A=G„corM/n 

A  =  D66  +  D22  +  GyzoorM2  1 "  ' 


A  =GyzcotM/n 

■^6  —  ( Gx/  core  ^yzcore)^ 


(3.6) 


Finite  Element  Modeling 


The  enhanced  version  of  SHELL  was  used  to  obtain  linear  and  nonlinear 
displacement  solutions  for  the  aforementioned  sandwich  structure.  Five  plates  with  the 
same  thickness  and  variable  widths  were  considered  in  order  to  compare  the  static 
responses  for  thin  and  thick  plates.  A  constant  value  of  h=12.192  mm  was  assumed,  and 
the  aspect  ratios  (S=a/h)  for  the  five  plates  were  10,  20,  30,  40  and  50.  The  meshes  for 
each  quarter-plate  model  consisted  of  N-by-N  square  elements  of  equal  size  (Ax=Ay= 
a  /2N).  The  appropriate  mesh  resolution  for  each  plate  was  obtained  through  convergence 
studies  of  meshes  ranging  from  12-by-12  to  24-by-24.  Furthermore,  the  maximum 
pressure  applied  to  each  plate  was  q0=6985  kPa. 

SHELL  cannot  automatically  consider  the  sinusoidal  load  distribution  in  Equation 
3.3.  Therefore  it  was  necessary  to  calculate  individual  nodal  loads.  This  was 
accomplished  by  integrating  the  product  of  Equation  3.3  and  the  Hermitian  shape 
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function  HK  ]  (associated  with  the  nodal  displacement  w  in  Equation  2.5  and  converted 
into  global  coordinates  based  on  Figure  2.3)  across  each  element’s  area  for  every  local 
element  node  K.  The  nodal  loads  for  an  entire  mesh  were  calculated  by  adding  the 
contributions  of  each  adjacent  element  for  a  given  node.  Note  that  this  method  ignores 
coupled  loads  obtained  from  the  other  shape  functions  in  Equation  2.5,  but  it  should  be  a 
good  approximation  of  the  distributed  load  for  meshes  with  a  relatively  large  number  of 
elements.  MATLAB™  [13]  and  Mathematica™  [22]  (commercial  mathematical  software 
packages)  were  employed  to  perform  the  integrations  and  other  computations.  The 
accuracy  of  its  nodal  load  calculations  was  verified  by  checking  for  load  symmetry  about 
each  quarter-plate’s  plane  of  symmetry  (x=y)  and  by  comparing  the  sum  of  all  nodal 
loads  to  the  total  force  from  a  continuous  distribution  (integrating  Equation  3.3  across  the 
entire  area  of  the  quarter-plate  yields  a  total  load  of  q0a  hi ). 


Convergence  Study 

As  shown  in  Table  3.1,  12x12  meshes  converged  well  for  S=10  and  20  while 
20x20  meshes  were  satisfactory  for  S=30,  40  and  50.  Deviations  under  5%  are 
considered  good  for  the  plate  elements  used  by  SHELL. 


Table  3.1 :  Displacement  Convergence  for  Plate  Meshes  (NL  Solution) 


wc  [mm]  for  NxN  mesh  at  q0=6985  kPa 

s 

12x12 

16x16 

20x20 

24x24 

%  Deviation 

10 

2.5038 

2.5058 

wmm  ■ 

SB 

20.853 

20.933 

HhH 

E9 

61.534 

61.783 

40 

111.86 

112.82 

0.858 

50 

167.29 

169.53 

1.34 

-4 
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Results  and  Discussion 


Figures  3.2  and  3.3  plot  the  nonlinear  FE  deflection  of  each  plate  at  its  center  (wc) 
due  to  varying  peak  load  pressures.  The  two  thickest  plates,  S=10  and  20,  exhibit  very 
linear  behavior  up  to  q0=6985  kPa.  On  the  other  hand,  geometric  nonlinearities  become 
evident  in  the  curves  for  S=30,  40  and  50  when  q0  is  greater  than  500  kPa.  At  the  highest 
load,  the  value  of  wc  for  the  thinnest  plate  (S=50)  is  about  15  times  greater  than  its 
thickness.  By  comparison,  a  linear  FE  solution  (from  Table  3.2)  for  S=50  would  predict  a 
maximum  deflection  that  is  58  times  greater  than  the  plate  thickness,  which  obviously 
violates  the  assumption  of  small  displacements  for  linear  behavior. 


Table  3.2:  Linear  FE  Results 


s 

q0  [kPa] 

wc  [mm] 

slope=  wc  /  q0  [mm/kPa] 

wc  [mm]  at  q0=6985  kPa 

10 

698.5 

0.251 

3.593e-4 

2.510 

m 

698.5 

2.313 

3.31  le-3 

23.13 

m 

698.5 

10.04 

1.437e-2 

100.4 

40 

698.5 

29.84 

4.272e-2 

298.4 

50 

698.5 

70.70 

1.012e-l 

707.0 

In  order  to  obtain  deflection  results  from  CLPT  and  Whitney’s  solution 
(Equations  3.5  and  3.6),  the  bending  stiffness  matrix  must  be  determined  for  each 
sandwich  plate.  SHELL’S  preprocessor  calculates  these  numbers  and  can  be  commanded 
to  display  them  in  the  output  file.  For  these  plates,  the  matrix  does  not  change  since  face 
and  core  thickness  are  held  constant.  The  relevant  values  (in  N-m=10  kPa-mm  )  are: 
Dn=13215.7 ,  D22=55 1.238  ,  D12=137.810  and  D66=272.014.  In  addition,  both  Pagano 
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and  Whitney  nondimensionalized  their  displacement  results  with  the  following 


expression: 


w  = 


100wcETTface 

q0S4h 


(3.7) 


Table  3.3  lists  the  nondimensional  displacements  from  CLPT,  Pagano’s  elasticity 
solution,  Whitney’s  solution  and  both  linear  and  nonlinear  FE  methods.  Figures  3.4,  3.5 
and  3.6  plot  these  values  versus  the  corresponding  aspect  ratio.  The  linear  FE  solution 
produced  deflections  similar  to  both  Pagano  and  Whitney  (albeit  slightly  stiffer  for  thick 
plates).  Furthermore,  all  linear  cases  converged  to  the  CLPT  solution  as  plates  became 
thinner.  The  differences  for  thick  plates  can  be  attributed  to  variations  in  how  each  theory 
considers  transverse  shear  effects  (or  neglects  them  in  the  case  of  CLPT).  As  bending 
effects  become  dominant  for  large  aspect  ratios,  each  linear  theory  produces  nearly 
identical  results.  The  divergence  between  linear  and  nonlinear  FE  results  as  S  is 
increased  is  due  to  a  coupling  between  bending  and  membrane  stiffness  that  is  not  present 
in  the  linear  case. 


Table  3.3:  Nondimensional  Plate  Deflection 


w 

s 

CLPT 

Pagano 

Whitney 

Linear  FE 

NL  FE 

q0=698.5  kPa 

NL  FE 

q„=3492.5  kPa 

NL  FE 

q0=6985  kPa 

1 

0.9238 

0.9238 

0.9238 

0.9238 

0.9238 

2.150 

1.300 

1.050 

0.950 

0.925 

2.535 

1.317 

1.076 

0.990 

0.950 

2.060 

1.186 

1.016 

0.956 

0.928 

2.060 

1.184 

0.996 

0.846 

0.658 

2.058 

1.148 

0.792 

0.501 

0.320 

2.054 

1.069 

0.623 

0.358 

0.220 
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Figure  3.6:  Nondimensional  Plate  Deflection  (Linear  and  NL)  vs.  Aspect  Ratio 


IV  Sandwich  Plate  Versus  Composite  Plate 

Finite  Element  Modeling 

In  order  to  demonstrate  some  of  the  advantages  of  sandwich  plates  in  terms  of 
stiffness-to-weight  ratios,  nonlinear  displacement  solutions  for  thin  and  thick  sandwiches 
were  obtained  from  SHELL  and  compared  to  solutions  (also  using  SHELL)  by  Owens 
[16,17]  for  similar  composite  plates.  His  six  plates  were  [02/90]s  laminates  of  a 
transversely-isotropic  graphite-epoxy  composite.  Each  ply  was  1 .016  mm  thick,  thus 
h=6.096  mm.  All  plates  were  square  and  their  widths  were  varied  to  produce  aspect 
ratios  (S=a/h)  of  10,  20,  30,  40,  50  and  60.  He  computed  nonlinear  finite  element 
solutions  for  simply-supported  edges  and  uniform  transverse  pressures. 

In  this  research,  each  sandwich  plate  had  the  same  overall  geometry  as  the  plates 
modeled  by  Owens  and  were  subjected  to  the  same  loads  and  boundary  conditions.  Each 
face  sheet  was  made  of  the  same  graphite-epoxy  material  in  a  [0/90 1/2]  lay-up  so  that  both 
faces  comprised  half  of  a  plate’s  total  thickness.  The  other  (central)  half  was  a 
honeycomb  core  made  of  Nomex™  (specifically  classified  as  HRH-10-1/8-9.0)  [8]. 
Figure  4.1  compares  the  geometry  of  both  types  of  plates.  The  walls  of  the  core  s 
hexagon-shaped  cells  run  parallel  to  the  Z-axis,  and  the  voids  between  the  cells  give  it 
negligible  in-plane  stiffness  compared  to  the  facesheets.  The  relevant  elastic  and  density 
properties  of  the  face  and  core  materials  are. 


Face:  Ell=1  37.9  GPa 

Glt=Glz=1.724  GPa 
vLt=0.25 
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Ett=3.447  GPa 
GTZ=0.6895  GPa 
pf =1.6  g/cm3  (from  [1]) 


(4.1) 


(4.2) 


Core:  ELL— Exx— GLT— 0 

GLZ=  120.6  MPa  Gxz=75.84  MPa 

vlx=0.5  pc  =0.14417  g/cm3 


Figure  4.1:  Composite  and  Sandwich  Quarter-Plate  Geometry 


Just  like  the  plates  in  the  Pagano /  Whitney  case  study,  the  quarter-plate  FE 
models  for  these  sandwiches  used  an  N  by  N  mesh  of  square  elements,  and  each  mesh 
was  refined  to  establish  convergence.  Furthermore,  each  plate  geometry  was  solved 
twice  to  obtain  results  for  core  orientations  of  [0°]  and  [90°].  This  was  done  to  see  if 
aligning  the  core’s  stiffest  transverse  plane  with  either  the  longitudinal  or  lateral 
directions  of  the  outer  face  plies  caused  significant  differences.  Finally,  the  maximum 
uniform  load  applied  to  each  plate  was  6985  kPa,  for  which  SHELL  automatically 
generated  the  equivalent  nodal  forces. 
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Convergence  Study 


Displacement  convergence  at  the  plate  center  for  the  highest  applied  pressure  was 
tested  for  three  of  the  six  plate  geometries.  The  results  are  listed  in  Table  4.1.  It  was 
found  that  12x12  meshes  were  acceptable  for  S=20  (and  S=10  since  it  has  a  stiffer 
response).  Similarly,  20x20  meshes  converged  reasonably  well  for  S=30,  40,  50  and  60. 


Table  4.1:  Displacement  Convergence  for  Sandwich  Plate  Meshes-[0°]  Core 


wc  [mm]  for  NxN  mesh  at  maximum  q0 

s 

12x12 

16x16 

20x20 

24x24 

Min.  %  Deviation 

12.327 

12.437 

0.892 

40 

48.417 

49.987 

51.778 

1.45 

60 

90.584 

94.399 

97.158 

99.268 

2.17 

Results  and  Discussion 

Figures  4.2  through  4.7  display  the  deflection  of  each  sandwich  plate’s  center  as  a 
function  of  applied  pressure,  and  Figures  4.8  and  4.9  combine  the  first  six  graphs  into  a 
family  of  curves  for  easier  comparison.  The  thickest  plates,  S=10,  behaved  very  linearly 
over  the  entire  pressure  range,  and  the  S=20  plates  were  linear  up  to  about  q0=2000  kPa 
and  then  showed  a  very  shallow  nonlinear  deviation.  On  the  other  hand,  the  thinner 
plates’  deflections  became  highly  nonlinear  at  pressures  less  than  1000  kPa.  In  addition,  a 
[90°]  orientation  of  the  core  caused  the  plate  to  be  slightly  more  flexible  for  S=10  and  20, 
but  for  thinner  plates,  the  curves  were  almost  identical.  In  fact,  Figures  4.5,  4.6  and  4.7 
show  only  the  results  for  a  [0°]  core  to  avoid  redundancy.  For  a  square  plate  with  the 
same  kind  of  support  on  every  side,  rotating  the  core  through  a  right  angle  has  the  same 
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effect  as  rotating  both  faces  instead  of  the  core.  It  is  important  to  keep  in  mind  that 
modifying  the  orientation  of  the  same  core  or  faces  should  have  very  different 
consequences  if  the  whole  plate  was  rectangular  or  had  multiple  types  of  edge  supports. 

Owens  published  deflection  results  [17]  in  a  nondimensionalized  form  for  selected 
pressures  (which  are  also  nondimensionalized)  as  a  function  of  aspect  ratio.  The 
nondimensional  forms  of  plate  center  deflection  and  applied  pressure  are: 

-  _  10wcELLface 

W  — - : - 

q0S4h  (4.3) 

q  =  10  q0  I  Ell  face 

While  q  is  directly  proportional  to  q0 ,  the  S4h  term  in  w  serves  to  cancel-out  geometric 
effects  on  wc  for  CLPT  solutions.  Deflection  results  from  the  previous  chapter  (Figures 
3.4  through  3.6)  showed  that  w  was  a  horizontal  line-  unaffected  by  changes  in  S—  for 
the  CLPT  case.  Furthermore,  the  weight  of  each  plate  can  be  considered  by  multiplying 
w  by  a  ratio  of  the  plate’s  overall  density  to  the  density  of  the  face  material  pf .  The 
composite  plates  obviously  had  densities  equal  to  pf ,  so  its  nondimensional 
displacements  were  unchanged.  On  the  other  hand,  the  core  accounts  for  half  of  each 
sandwich  plate’s  volume.  Therefore,  the  sandwich  plates  had  an  overall  density  equal  to 
the  average  of  the  core  and  face  densities.  From  the  data  in  Equations  4.1  and  4.2: 

Ps  =i(Pf  +  pc)  =  0.872  g/  cm3  =  0.545pf  (4.4) 

The  resulting  nondimensional  variable  wp  =  w  (ppiate  /  Pf)  can  be  used  as  an  indicator 
for  comparing  two  plates’  stiffness-to-weight  ratios. 
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Figures  4.10  through  4.17  compare  the  nondimensional  nonlinear  displacements 
of  each  composite  plate  and  its  corresponding  sandwich  plate  (the  same  value  of  S)  at 
four  nondimensional  load  levels.  Each  consecutive  pair  of  figures  plots  w  and  wp 
versus  S  for  a  given  q  .  In  Figures  4.10,  4.12,  4.14  and  4.16,  the  sandwich  constructions 
were  up  to  85%  more  flexible  than  the  composites  .  This  is  because  the  core  material 
provides  little  bending  stiffness  (none  for  this  FE  modeling)  and  offers  less  resistance  to 
transverse  shear  than  the  composite’s  additional  face  material  in  the  central  plies.  Hence, 
the  outer  faces  of  each  sandwich  must  compensate  by  bending  more.  Also  note  that  thick 
sandwiches  gain  flexibility  at  a  greater  rate  than  the  composite  plates  as  S  is  reduced. 

On  the  other  hand,  stiffness  contributions  due  to  bending  and  transverse  shear 
from  plies  near  the  midsurface  (either  face  or  core  material)  have  less  significance  in 
thinner  plates.  This  caused  the  sandwich  and  composite  curves  to  converge  sharply  at 
first  and  then  become  nearly  parallel  as  S  increased  and  bending  in  the  outer  faces 
became  dominant.  In  addition,  higher  loads  increased  the  rate  at  which  the  sandwich  and 
composite  curves  converge  and  decreased  the  ultimate  (large  S)  curve  deviation. 

When  the  nondimensional  deflections  of  the  sandwich  plates  were  scaled-down  to 
consider  their  lighter  weights,  the  resulting  plots  (Figures  4.1 1, 4.13, 4.15  and  4.17) 
suggested  equal  or  better  stiffness-to-weight  ratios  than  those  of  the  composite  plates  for 
a  wide  range  of  aspect  ratios.  The  wp  versus  S  curves  also  indicate  that  the  range  of  S  in 
which  the  particular  sandwich  outperforms  the  composite  is  bounded  due  to  merging  or 
crossing  of  the  curves.  Near  S=10,  the  curves  intersected  due  to  the  steeper  increase  in 
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flexibility  that  is  present  in  thick  sandwiches.  Furthermore,  the  sandwiches  have  lower 
axial  stiffness  for  response  to  nonlinear  membrane  and  flexural  coupling,  which  may 
cause  the  curves  to  intersect  again  beyond  S=60.  The  apparent  advantages  demonstrated 
by  these  sandwich  plates,  in  terms  of  specific  stiffness,  may  not  hold  true  for  all  sandwich 
and  composite  plate  constructions,  but  the  specific  case  illustrates  the  potential  benefits  of 
using  sandwich  plates. 


4-6 


2 


Figure  4.2:  Plate  Center  Deflection  vs.  Uniform  Pressure 
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Figure  4.3:  Plate  Center  Deflection  vs.  Uniform  Pressure 
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Figure  4.8:  Plate  Center  Deflection  vs.  Uniform  Pressure 


Figure  4.9:  Plate  Center  Deflection  vs.  Uniform  Pressure 
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Figure  4.1 1 :  Nondimensional  Plate  Center  Deflection  (incl.  weight)  vs.  Aspect  Ratio 
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Figure  4.14:  Nondimensional  Plate  Center  Deflection  vs.  Aspect  Ratio 


Figure  4.15:  Nondimensional  Plate  Center  Deflection  (incl.  weight)  vs.  Aspect 
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Figure  4.16:  Nondimensional  Plate  Center  Deflection  vs.  Aspect  Ratio 


Figure  4.17:  Nondimensional  Plate  Center  Deflection  (incl.  weight)  vs.  Aspect  Ratio 
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V.  Sandwich  Plate  Incipient  Damage  Predictions 


Finite  Element  Modeling 

The  preceding  displacement  solutions  for  various  sandwich  plates  were  obtained 
with  SHELL  under  the  assumption  of  perfect,  linear  elastic  materials.  However,  an  actual 
sandwich  may  experience  some  kind  of  initial  failure  well  before  any  geometric 
nonlinearities  take  effect.  Once  this  occurs,  any  solution  SHELL  generates  beyond  that 
point  is  invalid  due  to  the  presence  of  physical  nonlinearities  which  the  present  code 
cannot  consider.  Therefore,  FAILURE  was  written  for  the  purpose  of  attempting  to 
predict  where,  when  and  how  a  plate  initially  fails  (using  the  maximum  stress  criteria 
presented  in  Chapter  2).  Sample  plate  models  were  tested  to  verify  the  accuracy  of  its 
numerical  and  logical  procedures  (compared  to  manual  calculations  and  expected  output). 
Hence,  the  next  step  was  to  check  the  validity  of  its  methodology  by  modeling  actual 
plates  and  comparing  FAILURE’S  results  to  what  really  happens.  As  previously 
mentioned  in  Chapter  1 ,  a  low- velocity  impact  for  a  composite  plate  has  been  found  to 
have  internal  failure  characteristics  that  can  be  predicted  with  a  quasi-static  response 
[10,12,14],  Thus,  experimental  impact  studies  on  sandwich  plates  by  Harrington  [7] 
provided  a  means  for  applying  the  FAILURE  program. 

In  Harrington’s  work,  simply-supported  square  plates  (a=127  mm)  were  subjected 
to  impact  loads  at  their  center  by  the  dropping  of  a  spherical-nose  punch  from  a  series  of 
heights.  Four  cases  of  plate  and  load  combinations  were  considered  for  finite-element 
modeling,  and  these  are  listed  in  Table  5.1.  The  4  and  16-ply  facesheets  were  made  of 
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AS4/3501-6  graphite-epoxy  lamina  with  a  thickness  of  0.127  mm  per  ply.  The  4-ply 
faces  had  [0/90]s  lay-ups,  while  the  16-ply  faces  had  [0/90]4s  arrangements.  Each  plate’s 
core  was  12.7  mm  thick  and  made  of  an  HRX- 10-1/8-9.0  Nomex™  honeycomb  material 
[8].  Furthermore,  two  0.254  mm  thick  adhesive  layers  of  epoxy  bonded  each  face  to  the 
core.  For  convenient  referencing,  the  plates  with  4  and  16-ply  faces  are  denoted  as 
Sandwich  A  and  B,  respectively.  In  addition,  the  numbering  of  the  cases,  one  through 
four,  in  Table  5.1  corresponds  to  the  radii  of  plate  indentation  (R)  from  smallest  to 
largest.  For  each  case,  the  peak  load  (Pp)  equals  the  maximum  instantaneous  force  acting 
on  the  plate  within  the  experimental  impact  time  history,  and  it  is  assumed  quasi-static  for 
FE  analysis. 


Table  5.1 :  Plate  and  Loading  Cases 


Case 

Sandwich  Type 
(#  Face  Plies) 

Plate  Indentation  Radius 
from  Impact:  R  [mm] 

Peak  Load 
Pp[N] 

1 

B  (16) 

3.81 

3304.6 

2 

A  (4) 

5.08 

1620.9 

3 

B  (16) 

6.35 

3914.0 

4 

A  (4) 

7.62 

1969.7 

The  elastic  properties  of  the  face,  core  and  adhesive  layers  have  the  following 
relevant  values: 

Face:  ELL=1 19.3  GPa  ETT=9.098  GPa 

Glt=GI7=5.398  GPa  GTZ=4.319GPa  (5.1) 

vLX=0.25 

Core:  ELL=ETT=GLT=0 

Glz=120.6  MPa  GTZ=75.84  MPa  (5.2) 

vlt=0.5 
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Adhesive  (assumed  isotropic):  E=3.447  GPa 

v=0.35  (5.3) 

In  addition,  the  following  maximum  stress  values  were  known  or  assumed: 

Face:  aLLmax=2.016  GPa  aLLmin  = -1.398  GPa 

^TT max  =56.96  MPa  aTTmin  = -246.7  MPa  (5.4) 

^LT  max  ~®LZ  max  — 177.9  MPa 
°TZ  max  =142.3  MPa 

Core:  aZZmin  = -14.55  MPa  (5.5) 

&lz max  =177.9  MPa  aTZmax  =142.3  MPa 

Adhesive:  amax  uniaxial  =108.8  MPa  (5.6) 

Face  Ply  Interface:  axz  max  =aYZ  max  =142.3  MPa  (5.7) 

FAILURE  is  not  presently  designed  to  directly  use  the  core’s  transverse  compressive 
strength  (azz  min  from  Equation  5.5)  as  a  failure  condition.  This  is  because  the  user- 
defined  Z-coordinate  at  which  azz  is  calculated  may  or  may  not  be  the  location  of 
maximum  compressive  gzz  within  the  core.  However,  the  reported  czz  values,  if 
appropriate,  can  be  manually  used  to  check  for  initial  failure  due  to  crushing  of  the  core. 

The  impact  punch  was  assumed  to  create  an  ellipsoidal  pressure  distribution  on 
each  plate’s  top  surface  within  its  radius  of  indentation.  This  type  of  loading  is  similar  to 
that  obtained  from  Hertz  contact  stresses  for  isotropic  materials  [9].  Wu  and  Yen  [23] 
also  observed  the  presence  of  ellipsoidal  distributions  for  composite  plates  in  contact  with 
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rigid  spheres.  The  ellipsoidal  pressure  (q)  as  a  function  of  impact  force  (P)  and  radial 
distance  from  the  center  of  the  plate  (r),  has  the  following  expression: 


q(r)  =  q0(\-r2/R2)m 

r(x,y)  =  (x2+y2)m 
P=J  J  q(r)rdr  dQ  =  §7t  q0  R2 


Note  that  the  pressure  is  a  peak  value  of  q0  at  the  center  of  the  plate  (x  =  y  =  r  =  0)  and 
goes  to  zero  along  the  indentation  radius  (r  =  R).  Also  note  that  P  refers  to  the  impact 
force  on  the  entire  plate,  while  q0  has  the  same  pressure  for  both  full  and  quarter-plate 
models.  The  load  distribution  parameters  for  each  FE  case  are  listed  in  Table  5.2  when 
P=P 

r  i  p  . 


Table  5.2:  Ellipsoidal  Impact  Pressures 


Case 

R  [mm] 

Pp[N] 

Maximum  pressure  at  plate  center:  q0  [MPa] 

1 

3.81 

3304.6 

108.7 

2 

5.08 

1620.9 

30.00 

3 

6.35 

3914.0 

46.35 

4 

7.62 

1969.7 

16.20 

These  pressure  distributions  must  be  converted  to  discrete  nodal  loads  in  order  to 
use  SHELL.  The  method  employed  in  Chapter  3  (pp.  3-3  and  3-4)  was  also  used  here  by 
substituting  Equation  5.8  for  Equation  3.3  when  integrating  the  product  of  the  pressure 
and  shape  functions.  Figure  5.1  shows  three  mesh  arrangements  that  were  used  in  the 
elliptical  pressure  zones.  All  plate  centers  were  located  at  the  node  in  the  lower-left 
comer  of  each  mesh.  It  was  necessary  to  approximate  the  circular  indentation  zone  with 
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rectangles  because  FAILURE  is  limited  to  rectangular  meshes.  As  a  consequence,  those 
elements  which  lie  along  the  arc  had  to  be  integrated  over  a  partial  area  when  calculating 
their  nodal  loads.  In  order  to  minimize  the  number  of  elements  requiring  partial 
integration,  the  mesh  lines  in  Figure  5.1  were  designed  so  they  intersected  the  arc  at 
nodes.  The  3x3  mesh  was  employed  in  each  FE  case,  and  the  other  two  were  only  used  in 
case  4  as  part  of  the  convergence  study. 


3x3  zone 


2x2  zone 


lxl  zone 


Figure  5.1:  Quarter-Plate  FE  Meshes  in  Impact  Zone  (L;j  are  element  labels) 


Figures  5.2  through  5.9  display  the  entire  mesh  for  each  quarter  plate  model. 
Again,  the  center  of  each  plate  is  the  lower-left  corner  node.  The  size  of  each  element 
can  be  determined  from  the  element  increment  lengths  (the  distances  between  the  nodes) 
listed  in  Table  5.3-  which  are  the  same  for  both  the  X  and  Y-directions  (defined  upwards 
and  to  the  right  respectively)  since  all  plates  are  square.  The  four  24-by-24  meshes  were 
the  primary  models  used  for  checking  each  case  for  failure,  and  the  others  applied  to  case 
four  for  convergence  studies.  The  23-by-23  and  22-by-22  meshes  were  coarser  inside  the 
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impact  zone,  while  the  1 8-by- 1 8  and  12-by-12  meshes  were  coarser  outside  the  impact 


zone. 


Table  5.3:  Element  Sizing  of  Various  Meshes 


Sandwich/  Loading 
Case 

1 

2 

3 

4 

4 

4 

4 

4 

Mesh  Resolution 

24x24 

24x24 

24x24 

24x24 

23x23 

22x22 

18x18 

12x12 

Corresponding 

5.2 

5.3 

5.4 

5.5 

5.6 

5.7 

5.8 

5.9 

Figure 

Element  Increment 

Increment  Length  (N 

ode  Spacing)  | 

mm] 

1st 

1.905 

2.540 

3.175 

3.810 

5.388 

7.620 

3.810 

3.810 

2nd 

1.395 

1.859 

2.324 

2.789 

2.232 

0.707 

2.789 

2.789 

3rd 

0.511 

0.681 

0.851 

1.021 

0.707 

1.415 

1.021 

1.021 

4th 

0.756 

0.740 

0.723 

0.707 

1.415 

2.830 

0.707 

2.122 

5th 

1.511 

1.479 

1.447 

1.415 

2.830 

2.830 

1.415 

5.659 

6th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

8.489 

7th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

8.489 

8th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

5.659 

8.489 

9th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

5.659 

8.489 

10th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

5.659 

5.659 

11th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

5.659 

5.659 

12th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

5.659 

2.830 

13  th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

5.659 

14th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

15th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

16th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

17th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

18th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

2.830 

19th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

20th 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

21st 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

22nd 

3.023 

2.958 

2.894 

2.830 

2.830 

2.830 

23rd 

3.023 

2.958 

2.894 

2.830 

2.830 

24th 

3.023 

2.958 

2.894 

2.830 

Note:  1st  Increment  is  at  the  center  of  the  plate  (lower  left  comer  of  each  mesh) 
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Figure  5.2:  Case  1,  24x24  Mesh 


Figure  5.3:  Case  2,  24x24  Mesh 
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Figure  5.5:  Case  4, 24x24  Mesh 
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Figure  5.6:  Case  4,  23x23  Mesh 


Figure  5.7:  Case  4, 22x22  Mesh 
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Figure  5.9:  Case  4, 12x12  Mesh 


5-10 


Convergence  Study 


The  24-by-24  mesh  for  case  4  (Figure  5.5)  underwent  two  types  of  refinement  to 
check  for  convergence  of  the  displacement  results  (within  5%  relative  deviation).  First, 
the  element  spacing  in  regions  outside  the  load  zone  was  widened  to  form  1 8-by-l  8  and 
12-by-12  meshes  (Figures  5.8  and  5.9).  The  profile  of  plate  deflection  along  the  line  x=y 
was  compared  for  each  mesh  at  the  peak  load.  The  profiles  turned  out  to  be  identical 
(very  small  numerical  deviations  could  not  be  identified  graphically).  This  suggests  that 
these  plate  models  were  practically  insensitive  to  changes  away  from  the  applied  load. 

The  profile  is  not  included  here,  since  it  does  not  show  anything  extraordinary  in 
comparison  to  a  typical  profile  for  a  simply-supported  plate  with  a  central  transverse  load. 

Second,  the  mesh  arrangement  within  the  load  zone  for  case  4  was  modified 
according  to  Figure  5.1.  This  produced  the  23-by-23  and  22-by-22  plate  meshes  shown  in 
Figures  5.6  and  5.7.  Calculations  for  azz  were  obtained  for  each  mesh  at  the  top  of  the 
plate  near  the  center  (where  the  pressure  is  highest).  The  distance  from  the  center  of  the 
plate  to  the  closest  Gauss  points  in  local  elements  Ln  ,L2\  and  L31  (Figure  5.1)  varies 
slightly  for  each  mesh  (Gauss  point  positions  are  fixed  in  the  natural  coordinate  system 
and  scaled  for  each  element’s  dimensions).  However,  they  are  all  near  enough  to  the 
center  that  the  variation  in  the  elliptical  pressure  distribution  is  negligible  with  respect  to 
the  peak  value.  For  q0=16.20  MPa,  the  calculated  values,  in  decreasing  order  of  mesh 
resolution,  were:  16.14 , 13.09  and  10.46  MPa.  Hence,  the  24-by-24  mesh  converged  to 
within  0.4%  of  the  actual  value. 
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Results  and  Discussion 


For  every  finite  element  case,  the  plate’s  maximum  deflection  at  peak  load  was 
very  small  relative  to  its  geometry.  Sandwich  A  (4  ply  face)  is  the  more  flexible  plate, 
and  for  it’s  highest  load  (case  4)  the  midsurface  deflection  at  the  plate’s  center  was  only 
0.71mm—  about  5%  of  its  total  thickness  (hA=  14.224  mm).  Therefore,  although  each 
case  was  solved  nonlinearly,  its  static  results  were  practically  linear. 

Table  5.4  lists  the  incipient  failure  results  for  each  case,  in  which  all  maximum 
stresses  except  transverse  core  compression  were  considered  (Equations  5.4  through  5.7). 


Table  5.4:  Initial  FE  Sandwich  Failure  Ignoring  Core  Compression 


Case  (Sandwich 

Type  /  #  Face  Plies) 

Pp[N] 

%  of  Pp  Range  Mode  of  Location  (material) 

for  First  Failure  Failure 

1  (B/16) 

3304.6 

110-120  Lateral  Bottom  of  plate  near  center 

Tension  (face) 

2  (A/4) 

1620.9 

80-90  Transverse  Midsurface  near  center  (core) 

Shear 

3  (B/16) 

3914.0 

150-160  Lateral  Bottom  of  plate  near  center 

Tension  (face) 

4  (A/4) 

1969.7 

100-110  Transverse  Midsurface  near  center  (core) 

Shear 

In  each  case,  the  initial  failure  occurred  near  or  beyond  the  peak  impact  force.  Tension 
fracture  of  the  face’s  matrix  material  on  the  bottom  of  the  plate  is  a  possible  failure  mode 
for  a  standard  laminated  composite,  but  it  is  not  realistic  for  these  kinds  of  sandwiches. 
The  lower  surface  of  the  top  face  is  more  likely  to  fail  that  way  as  a  progressive  mode 
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when  core  damage  from  indentation  removes  the  top  face’s  localized  support  in  the  Z- 
direction.  It  occurs  on  the  bottom  face  for  sandwich  B  because  SHELL’S  plate  kinematics 
(Equation  2.1)  prevents  describing  indentation.  The  top  and  bottom  faces  (for  a  fixed  x 
and  y)  must  have  the  same  w-translation  as  the  midsurface.  Transverse  shear  failure  of 
sandwich  A’s  core  is  a  possible  mode  in  later  stages,  but  core  crushing  is  more  likely  to 
occur  first. 

Before  examining  the  core  for  transverse  compression  failure,  it  was  necessary  to 
see  if  Equation  2.12  was  successful  in  obtaining  values  of  ctzz  in  the  impact  zone  that 
were  close  to  the  applied  pressure  (Equation  5.8)  on  the  top  surface  (or  at  least  close  to  q0 
near  the  plate’s  center).  Table  5.5  compares  the  known  pressures  to  the  calculated 
stresses. 


Table  5.5:  Maximum  Compressive  Stress  Results  at  Top  and  Center  of  Plate 


Case 

q0  [MPa] 

aj.0,0,  -h/2)  [MPa] 

%  error 

1 

108.7 

70.44 

35 

2 

30.00 

28.37 

5.4 

3 

46.35 

45.72 

1.3 

4 

16.20 

16.14 

0.4 

For  cases  2  through  4,  the  stress  results  were  very  good,  but  case  1  underpredicted  the 
applied  pressure  by  a  large  margin.  There  is  an  inverse  correlation  between  the  error  and 
radius  of  the  pressure  zone  (the  indentation  radius).  Since  pressure  must  drop  from  a  peak 
value  to  zero  within  this  zone,  its  gradient  may  be  too  high  in  case  1  for  the  employed 
mesh  .  However,  finer  meshes  were  not  generated  because  the  predominantly  manual 
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process  used  to  calculate  the  nodal  loads  limited  the  practicality  of  finer  meshes. 
Mathematica™  [22]  was  used  for  the  complex  integrations,  but  due  to  difficulties  in 
automating  the  entire  process,  each  element’s  nodal  loads  had  to  be  individually 
evaluated  and  manually  added  to  the  nodal  loads  in  adjacent  elements.  Besides,  the  other 
three  cases  provided  sufficient  data  for  comparing  core  crushing  predictions. 

Table  5.6  lists  each  case’s  calculated  value  of  azz(0,0,  -hc/2)  at  the  top  of  the  core 
and  the  center  of  the  plate,  using  Equation  2.12,  for  the  highest  load  levels.  The  linear 
behavior  exhibited  by  all  the  plates  permitted  the  use  of  linear  interpolation  to  determine 
the  failure  loads  at  which  the  FE  models  predicted  core  crushing. 

Table  5.6:  FE  Failure  Load  Predictions  for  Core  Compression 


Compared  to  the  other  initial  failure  results  in  Table  5.4,  crushing  of  the  core  within  the 
impact  region  occurred  at  substantially  lower  fractions  of  the  maximum  load  for  each 
case,  except  4.  Furthermore,  although  case  4  did  not  show  failure  from  any  mode  until 
the  peak  load  was  exceeded,  both  tables  predicted  some  kind  of  core  failure  (either 
crushing  or  shearing)  at  nearly  the  same  load  level.  Hence,  through  the  use  of  an 
analytical  method  that  assumed  zero  transverse  normal  stresses  in  determining  the  static 
response,  it  was  still  possible  to  detect  transverse  compression  of  the  core  as  the  primary 


Case 

q0  [MPa] 

1 

108.7 

2 

30.00 

3 

46.35 

^(0,0,  -hc/2) 

[MPa]  at  q0 

%  of  q0  for  core  failure  at  14.55  [MPa] 

pct=14.55  /  ctzz(0,0,  -hc/2)  x  100% 

Pp[N] 

FE  Failure 

Load  [N] 

Pp  x  pet 

42.94 

34% 

3304.6 

1123.6 

23.56 

62% 

1620.9 

1005.4 

26.26 

55% 

3914.0 

2152.7 

13.40 

109  % 

1969.7 

2146.9 

mode  (or  one  of  several  modes)  in  which  these  types  of  sandwich  plates  initially  fail. 
Regardless  of  the  projected  load  levels,  the  mode  agrees  with  Harrington’s  experimental 
findings  which  showed  incipient  indentation  damage,  at  least  partially  due  to  core 
crushing,  within  the  impact  zone  for  sandwiches  with  4  and  16-ply  faces. 

In  Figures  5.10  through  5.13,  the  FE  load  levels  at  predicted  core  failure  have 
been  superimposed,  for  their  respective  cases,  onto  the  actual  time  histories  of  impact 
loading  from  Harrington’s  experiments.  Actual  initial  failure  is  usually  represented  by 
the  first  sharp  drop  in  load  (in  excess  of  noise  on  the  curve),  which  signifies  a  sudden 
shift  in  a  plate’s  equilibrium  state  due  to  a  reduction  in  stiffness.  For  the  16-ply  face 
sandwiches  (type  B)  in  cases  1  and  3,  such  drops  were  clearly  evident  at  loads  of  about 
2500  N  for  both  Figure  2. 10  and  2. 12.  The  FE  failure  load  in  case  1  underestimated  the 
onset  of  failure  by  about  50%,  but  case  3  provided  a  very  close  estimate  of  the  actual 
failure  load. 

For  the  4-ply  face  sandwiches  (type  A)  in  cases  2  and  4,  the  actual  failure  loads 
were  not  clearly  distinguishable  from  the  noise  present  in  Figures  5.11  and  5.13,  although 
moderate  spikes  between  500  and  1000  N  may  or  may  not  be  due  to  failure.  The 
sandwich  A  plates  experienced  more  data  noise  than  the  sandwich  B  plates  because  they 
were  more  flexible  and  thus  subjected  to  greater  momentum  transfer  during  the  impact.  In 
addition,  sandwich  A  displayed  the  same  phenomenon  as  sandwich  B  in  the  finite 
element  results-  a  near  doubling  of  the  projected  first-core-crushing  load  for  the  case 
with  the  higher  applied  load.  One  would  expect  the  loads  that  cause  initial  internal  failure 
to  remain  constant  for  different  applied  loads  on  the  same  plate  (as  they  clearly  do  for 
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sandwich  B).  This  inconsistency  between  the  finite  element  and  experimental  results 
suggests  that  maximum  stress  may  not  be  the  best  choice  of  failure  criteria.  Although  the 
primary  mode  was  core  crushing,  other  stress  components  may  contribute  to  the  initial 
failure.  The  maximum  stress  criteria  isolate  the  stresses  and  do  not  permit  coupling  to 
affect  the  failure  predictions. 

Another  source  of  inconsistency  in  the  failure  results  could  be  the  FE  modeling  of 
the  pressure  distribution.  Cases  2  and  4  had  larger  impact  forces  than  cases  1  and  3, 
respectively,  but  their  larger  indentation  radii  (the  assumed  constant  radii  of  the  impact 
zones)  spread  out  the  loading  so  that  the  peak  pressures  actually  dropped  for  higher  total 
loads,  as  shown  in  Table  5.2.  In  the  former  cases,  the  values  of  q0  were  roughly  cut  in 
half,  which  explains  the  doubling  of  the  projected  load  for  initial  internal  failure.  In  an 
actual  impact,  both  total  load  and  contact  area  vary  with  time,  and  the  right  combination 
produces  a  peak  pressure  high  enough  to  initiate  failure  in  the  core.  This  type  of 
nonlinear  behavior  cannot  be  modeled  with  SHELL;  therefore,  the  good  results  that 
occurred  in  case  3  were  most  likely  a  coincidence  brought  about  by  nearly  having  the 
right  mesh  size  to  produce  the  actual  failure  load. 
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VI.  Conclusions 


In  this  thesis  a  geometrically  nonlinear  finite  element  program,  created  for  static 
analysis  of  composite  plates  and  shells,  was  enhanced  so  that  it  could  be  used  to  study 
sandwich  plates.  Furthermore,  a  new,  separate  postprocessing  unit  was  created  in  order 
to  detect  initial  failure  in  a  plate  using  the  maximum  stress  criteria.  The  program  was 
also  given  the  capability  of  estimating  the  transverse  normal  stresses  within  a  plate  by 
enforcing  equilibrium  through  its  thickness.  Three  case  studies  were  investigated  for  the 
following  purposes: 

1.  To  validate  the  sandwich  plate  enhancements  to  the  FE  code  by  comparing  its 
displacement  results  to  those  of  established  linear  solutions  for  a  particular  sandwich 
plate  problem. 

2.  To  compare  the  stiffness  and  stiffness-to-weight  characteristics  of  regular 
composite  plates  to  those  of  sandwich  plates  for  different  load  intensities  and  aspect 
ratios. 

3.  To  simulate  low- velocity  impact  tests  on  sandwich  plates  with  a  quasi-static 
FE  solution  and  attempt  to  predict  incipient  plate  damage  using  the  maximum  stress 
criteria. 
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Linear  solutions  for  simply-supported  sandwich  plates,  under  a  sinusoidal 
transverse  pressure,  showed  good  agreement  with  Pagano’s  exact  elasticity  solution  and 
Whitney’s  laminated  plate  solution  for  both  thick  and  thin  plates.  For  thin  plates,  all  three 
methods  converged  to  the  CLPT  solution.  The  code  enhancements  related  to  sandwich 
plates  only  affected  the  formulation  of  the  constitutive  relations  in  the  preprocessor.  Since 
they  do  not  change  for  either  a  linear  or  nonlinear  solution  of  the  same  plate,  the  modified 
code  can  be  considered  valid  for  sandwich  plates  using  either  solution  method. 

Comparisons  between  a  graphite-epoxy  composite  and  a  sandwich  with  similar 
facesheets,  a  honeycomb  core  and  the  same  overall  geometry,  show  that  the  sandwich  is 
more  flexible  (especially  for  thick  plate).  The  differences  in  stiffness  become  smaller  for 
thin  plates  as  bending  in  the  outer  faces  dominate  the  response.  When  the  lighter  weight 
of  the  core  material  is  taken  into  account,  the  sandwich  plate  demonstrates  a  significantly 
higher  stiffness-to-weight  ratio  than  the  composite  for  both  thick  and  thin  plate  within 
certain  bounds.  If  specific  stiffness  is  the  primary  criterion  in  selecting  a  material,  a 
sandwich  construction  may  be  a  better  choice  than  a  laminated  composite  of  similar 
construction. 

The  procedure  for  calculating  crzz  was  shown  to  be  capable  of  highly  accurate 
estimates  of  transverse  pressures  on  the  top  surface  of  a  plate,  provided  the  FE  mesh  in 
these  areas  were  properly  refined.  Therefore,  the  method  was  partially  successful  in 
extracting  a  three-dimensional  stress  state  from  a  two-dimensional  solution.  In  addition, 
for  the  case  of  sandwich  plates  subjected  to  low- velocity  impact  loads,  the  use  of  quasi- 
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static  FE  modeling  and  the  maximum  stress  criteria  was  successful  in  detecting  core 
crushing  in  the  impact  zone  as  one  of  the  primary  modes  of  incipient  damage.  However, 
inconsistencies  in  the  predicted  load  levels  for  initial  failure  suggest  the  need  for  a  more 
complex  criteria  that  considers  stress  coupling.  Although  the  presence  of  time-dependent 
nonlinearities,  like  momentum  transfer  and  variable  contact  areas,  also  contributed  to 
preventing  the  quasi-static  models  from  making  good  predictions  of  when  initial  failure 
occurred,  this  quasi-static  approach  was  at  least  able  to  identify  where  and  how  it 
occurred. 
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Appendix  A: 

Modifications  to  Finite  Element  Code  SHELL 


The  original  version  of  SHELL  was  written  by  Dennis  [4],  It  is  a  FORTRAN  code 
with  the  following  components: 


Table  A.  1 :  Components  of  SHELL  code 


Component  Name 

Description 

Main  program  SHELL 

Subroutines: 

Preprocessor,  solver  and  postprocessor 

MESH 

Automatically  generates  rectangular  mesh 

ELAST 

Calculates  constitutive  relations  and  elasticity  matrices 

STIFF 

Calculates  element  stiffness  matrix  and  force  vector 

SHAPE 

Evaluates  shape  functions  at  Gauss  points 

DIS 

Evaluates  displacement  gradients  at  Gauss  points 

BNDY 

Imposes  prescribed  boundary  conditions 

SOLVE 

Solves  linear  equation  systems 

STRESS 

Calculates  stresses  at  outer  Gauss  points  from  displacement  results 

INCREMENT 

Increments  loads  and/or  displacements  for  nonlinear  analysis 

CONVERGE 

Tests  convergence  of  nonlinear  analysis 

PK 

Calculates  independent  element  stiffness  matrix  K  for  a  plate 

PN1 

N 1  for  a  plate 

PN2 

N2  for  a  plate 

SK 

K  for  a  shell 

SN1 

N 1  for  a  shell 

SN2 

N2  for  a  shell 

STOP 

Stops  program  if  unable  to  converge  itterations  using  Riks  method 

CHSIGN 

Used  in  Riks  method  to  allow  backwards  incrementing 

The  code  required  three  major  modifications  for  use  in  this  research: 

1 .  Enhance  the  preprocessor  to  allow  sandwich  constructions 

2.  Restore  load  and  displacement  control  options  to  the  nonlinear  solver 

3.  Generate  a  secondary  output  file  for  use  with  a  separate  initial-failure  check 

program 

Since  the  entire  SHELL  code  is  very  large,  only  those  components  requiring  extensive 
changes  are  listed  at  the  end  of  this  appendix.  The  other  components  remain  unaltered  or 
needed  minor  corrections  in  its  nonexecutable  statements  (i.e.  altering  common  variable 
blocks  to  make  them  consistent  with  the  rest  of  the  program).  The  new  structure  for 
SHELL’S  input  deck  is  also  included  just  before  the  code  listings. 

Sandwich  Construction 

The  previous  version  of  SHELL  could  model  an  isotropic  material  or  a  symmetric 
laminate  consisting  of  plies  of  the  same  orthotropic  material  at  different  orientations.  It 
also  required  each  ply  to  have  the  same  uniform  thickness.  The  input  deck  for  a  laminate 
contained:  the  number  of  plies,  the  ply  thickness,  the  orthotropic  elastic  properties  and 
each  ply’s  orientation.  In  order  to  model  sandwich  constructions,  the  code  was  altered  to 
allow  multiple  sets  of  elastic  properties  and  variable  ply-to-ply  thicknesses  (although  still 
uniform  across  a  single  ply).  The  new  input  deck  for  a  laminate  includes:  the  number  of 
plies,  the  number  of  materials,  an  indicator  for  uniform  or  variable  ply  thickness,  each 
material’s  orthotropic  elastic  properties,  and  each  ply’s  orientation,  material  reference 
number,  and  thickness  (or  a  single  thickness  value  if  uniform).  Note  that  a  sandwich 
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containing  isotropic  materials  can  be  modeled  as  a  laminate  by  treating  their  elastic 
properties  as  orthotropic  but  numerically  consistent  with  isotropic. 

The  use  of  multiple  elastic  property  sets  was  implemented  into  the  code  by 
converting  the  single-value  variables  for  ELL ,  ETT ,  GLT ,  GLZ  ,  GTZ  and  vLT  into  one¬ 
dimensional  arrays  .  Since  plies  can  be  at  different  orientations,  the  preprocessor  was 
already  designed  to  calculate  a  separate  set  of  constitutive  relations  for  each  one. 
Therefore,  all  that  was  needed  was  a  way  to  index  the  correct  element  in  each  elastic 
property  array  (corresponding  to  the  ply’s  material  reference  number).  This  was  done  by 
creating  a  material-stacking-sequence  (MSS)  array  similar  to  the  preexisting  orientation 
angle  array.  As  the  code  cycles  through  each  ply,  it  reads  a  new  material  number  and 
angle  and  uses  the  former  to  load  the  proper  elastic  constants. 

Enhancing  the  code  to  allow  variable  ply  thicknesses  was  not  essential,  but  it 
could  greatly  reduce  redundancies  in  the  input  deck  and  calculations  throughout  the 
program.  For  example,  a  3 -ply  laminate  with  ply  thicknesses  of  5,  36  and  5  units  would 
otherwise  require  a  46-ply  model  with  unit  thickness  per  ply.  The  large  variation  between 
face  and  core  thicknesses  in  typical  sandwiches  would  amplify  this  redundancy.  The  old 
preprocessor  used  the  uniform  ply  thickness  to  calculate  the  through-the-thickness  (Z) 
coordinates  of  each  ply  at  its  upper,  middle  and  lower  surfaces.  A  modified  method  using 
a  ply  thickness  array  (similar  in  form  to  the  MSS  and  orientation  arrays)  is  now  employed 
when  variable  thickness  is  indicated.  The  method  involves  simple  step-by-step  addition 
and  is  too  elementary  to  warrant  an  in-depth  explanation. 
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Some  core  materials  have  negligible  stiffness  in  the  in-plane  directions.  However, 
setting  Ell=Ett=0  in  the  input  deck  will  cause  the  program  to  crash.  The  error  is  due  to  a 
line  that  calculates  vTL  from  the  relation: 


To  prevent  division  by  zero,  a  precautionary  step  was  added  to  the  code  which  sets 
vtl=vlt  and  skips  Equation  A.  1  whenever  ELL=ETT. 


The  original  version  of  SHELL  included  both  load  and  displacement  control 
methods  for  the  nonlinear  solver.  A  later  version  included  the  Riks-Wemper  method 
which  allows  a  better  description  of  a  cylindrical  shell’s  behavior  when  it  undergoes 
snapping  instability  [18  ].  However,  this  research  uncovered  the  fact  that  in  adding  the 
Riks  method  the  other  methods  had  been  removed.  Load  or  displacement  control  is 
usually  more  practical  when  modeling  flat  plates  because  plates  do  not  tend  to  snap,  and 
the  Riks  method  normally  does  not  increment  the  loads  and/or  displacements  at  regular 
intervals. 


In  converting  the  program  to  the  Riks  method,  certain  parts  of  the  original  code 
were  deleted  or  turned  into  comments.  Fortunately,  another  finite  element  program  called 
ISHELL  contains  a  processor  unit  nearly  identical  in  format  to  SHELL  but  features  the 
original  load  and  displacement  control  instead  of  the  Riks  method.  The  process  of 
enhancing  SHELL  to  include  all  three  methods,  involved  a  line-by-line  comparison  and 
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merging  of  both  processor  codes.  An  indicator  was  added  to  the  input  deck  to  trigger  the 
use  or  disuse  of  the  Riks  method,  and  many  if-then  statements  were  added  to  the  code  to 
skip  unwanted  operations  in  either  case. 

Secondary  Output  File 

In  order  to  execute  the  separate  initial-failure  check  program  FAILURE,  certain 
model-dependent  information  is  needed  from  SHELL.  The  required  data  includes: 

1 .  Model  parameters  from  the  input  deck  (an  isotropic  or  laminate  model,  a 

linear  or  nonlinear  solution,  the  number  of  elements,  plies  and  materials 
and  the  elastic  material  properties) 

2.  Preprocessor  calculations  (the  nodal  coordinates,  nodal  connectivity  of 

elements,  constitutive  relations,  z-coordinates  and  thickness  factor  for 
transverse  shear) 

3.  Nodal  displacement  results  (for  each  increment  if  solution  is  nonlinear) 

All  of  the  above  is  written  to  a  separate  output  file  in  a  format  easily  read  by  FAILURE. 

A  new  indicator  in  the  SHELL  input  deck  allows  the  user  to  decide  whether  or  not  to 
generate  the  file. 

Other  Modifications 

In  addition  to  the  aforementioned  code  changes,  several  other  nonessential 
modifications  were  implemented  to  enhance  the  user-friendliness  of  the  software.  First, 
direct  keyboard  input  was  added  to  the  beginning  of  the  program  to  allow  user-defined 
names  for  the  input  and  output  files.  This  alleviates  the  task  of  renaming  or  relocating  old 
files  before  running  a  new  model  with  default  filenames.  It  also  allows  simultaneous 
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execution  of  multiple  models  on  a  computer  network  without  the  need  for  extra  copies  of 
the  program  in  separate  file  directories. 

Second,  the  double-precision  variables  that  contain  nodal  coordinate  and 
displacement  values  were  reformatted.  Each  value’s  scientific  notation  is  now  printed  to 
the  output  files  with  an  “E”  (instead  of  a  “D”)  to  indicate  the  exponent.  This  data  can  be 
cut-and-pasted  into  separate  files  for  use  with  commercial  math  or  graphing  software. 
However,  it  was  discovered  that  some  software  packages  do  not  recognize  “D”  as  an 
acceptable  substitute  for  “E”  and  will  misread  the  data  or  generate  a  syntax  error.  The 
modified  output  format  alleviated  this  problem. 
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Card  Variable  Type  Variable  Description  &  Allowable  Contents/Array  Size  Notes 
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Appendix  B: 


Initial  Failure  Postprocessor  Code  FAILURE 

FAILURE  is  a  FORTRAN  code  for  use  with  the  updated  version  of  SHELL.  It  is 
an  enhanced  postprocessor  designed  to  predict  initial  failure  of  flat  plates  using  maximum 
stress  criteria.  It  also  features  a  procedure  for  estimating  maximum  transverse  normal 
stress  magnitudes  through  a  plate’s  thickness.  It  is  limited  to  finite  element  plate  models 
constructed  from  28  degree-of-freedom  elements  in  a  regular  rectangular  mesh.  Execution 
of  FAILURE  requires  two  input  decks.  The  first  is  generated  by  SHELL  for  a  given 
model,  and  the  second  is  user-defined.  The  structure  of  the  second  deck  and  a  sample  of 
program  output  is  contained  in  this  appendix  (just  before  the  listing  of  the  FAILURE 
code). 

FAILURE  consists  of  the  main  program  and  four  subroutines:  FSTRESS, 
PARFIT,  SHAPE  and  DIS.  The  latter  two  are  identical  to  those  used  by  SHELL  (see 
Appendix  A,  Table  A.l).  In  addition,  FSTRESS  is  a  modified  version  of  SHELL’S  own 
STRESS  subroutine.  The  main  difference  between  the  two  is  that  FSTRESS  does  not 
place  stress  calculations  in  temporary  variables  for  immediate  reporting  to  the  output  file. 
Instead,  it  stores  the  values  for  a  single  element  and  solution  increment  in  arrays,  so  they 
can  be  mathematically  manipulated.  Finally,  PARFIT  is  a  set  of  linear  algebra 
computations  used  to  curve  fit  discrete  stress  and  Z-coordinate  data  to  parabolic 
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distributions.  This  is  needed  for  stress  averaging  through  a  ply’s  thickness  when  checking 
material  failure  and  for  evaluating  azz  by  enforcing  equilibrium. 

The  program  gives  the  user  the  option  of  checking  for  failure  in  a  single  element 
or  all  of  them.  The  same  option  applies  to  load  or  displacement  control  increments  (one 
or  all)  in  a  nonlinear  solution.  All  criteria  stresses  should  be  entered  as  absolute  values 
since  negative  stresses  are  accounted  for  when  the  program  checks  for  failure. 
Furthermore,  any  of  the  material  or  shear  delamination  failure  modes  can  be  turned  off  by 
entering  a  value  0.0  for  its  maximum  stress  value. 

When  the  algorithm  for  evaluating  ctzz  (denoted  cr33  in  some  parts  of  the  code)  is 
activated,  the  user  must  define  a  single  Z-coordinate  to  identify  the  plate  surface  (on  or 
between  -h/2  and  h/2)  at  which  values  are  to  be  calculated.  In  addition,  a  nonzero  shear- 
to-moment  ratio  may  be  entered  in  order  to  limit  azz  calculations  to  elements  in  which 
transverse  shear  stresses  are  significant  compared  to  in-plane  bending  stresses.  For  this 
particular  research,  shear-moment  ratios  were  not  employed  (although  they  were  initially 
proposed  which  led  to  their  inclusion  in  the  code). 
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User-defined  Input  Deck  to  FAILURE 
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Sample  Program  Output 


Qtr  sandwich  case2  24x24mesh  4  ply  face  r=.2  lo  load 
Sandwich  model  case2  4ply  lo  load  sig33  at  top  of  plate 

INITIAL  FAILURE  CHECK 

ELEMENT  TO  CHECK  (0  FOR  ALL):  1 
INCREMENT  TO  CHECK(0  FOR  ALL):  15 
REPORT  AVERAGED  AND  TRANSFORMED  STRESSES 

FOR  FAILED  REGIONS  (1  FOR  YES,0  FOR  NO):  1 

NUMBER  OF  MATERIALS=  3 
NUMBER  OF  PLIES=  11 

MATERIAL  FAILURE  MODES/CRITERIA 
ORTHOTROPIC  (MAX  STRESS) 

1-  LONGITUDINAL  (FIBER)  TENSION 

2- LONGITUDINAL  (FIBER)  COMPRESSION 

3- LATERAL  (MATRIX)  TENSION 

4- LATERAL  (MATRIX)  COMPRESSION 

5- LONG/LAT  IN-PLANE  SHEAR 

6- LAT/Z  TRANSVERSE  SHEAR 

7- LONG/Z  TRANSVERSE  SHEAR 
ISOTROPIC  (MAX  SHEAR  STRESS) 

1 - UNIAXIAL  TENSION 

2- UNIAXIAL  COMPRESSION 

3- SHEAR 

CRITERIA  VALUES  (MAGNITUDES) 

MATERIAL  #  1  (ORTHOTROPIC) 

MODE  1:0. 292393E+06 
MODE  2:0.202778E+06 
MODE  3:0. 8261 00E+04 
MODE  4:0.357770E+05 
MODE  5:0.258000E+05 
MODE  6:0.206400E+05 
MODE  7:0.258000E+05 
MATERIAL  #  2  (ORTHOTROPIC) 

MODE  1 :0.000000E+00 
MODE  2:O.OOOOOOE+O0 
MODE  3:0.000000E+00 
MODE  4:0.000000E+00 
MODE  5:0.000000E+00 
MODE  6:0.300000E+03 
MODE  7:0.515000E+03 
MATERIAL  #  3  (ISOTROPIC) 

MODE  1 :0.157800E+05 
MODE  2:0. 157800E+05 
MODE  3:0.789000E+04 

INTER-PLY  SHEAR  DELAMINATION 
X-DIRECTION  (13) :0.206400E+05 
Y-DIRECTION  (23) :0.206400E+05 

ESTIMATE  MAX  (MAGNITUDE)  TRANSVERSE  NORMAL 
STRESS  (1  FOR  YES,0  FOR  NO):  1 
MINIMUM  (ABS  VALUE)  S I GMA23/S I GMA22  OR 
SIGMA13/SIGMA1 1  RATIO  FOR  TRANSVERSE 
NORMAL  ESTIMATE=  0.000000E+00 
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INITIAL  FAILURE  RESULTS 


INCREMENT  #  15 

ELEMENT  #  1 

PLY  #  6  MATERIAL  #  2 

MATERIAL  FAILURE  MODES 
6 


AVERAGED  PLY  STRESSES  (X,Y,S1 1 ,S22,S12,S23,S13) 

.25446E-02  .25446E-02  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.950280E+01 

.25446E-02  .97455E-01  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.352938E+03 

.97455E-01  .25446E-02  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.873531E+01 

.97455E-01  .97455E-01  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.323595E+03 

TRANSFORMED  PLY  STRESSES  <X(Y,SLL,STT,SLT,STZ,SLZ) 

. 25446E-02  .25446E-02  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.950280E+01 

.25446E-02  .97455E-01  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.352938E+03 

. 97455E -01  .25446E-02  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.873531E+01 

.97455E-01  .97455E-01  O.OOOOOOE+OO  O.OOOOOOE+OO  O.OOOOOOE+OO  -.323595E+03 


- . 1 19206E+02 
- . 107059E+02 
-. 438665 E+03 
- .392255E+03 

- . 1 19206E+02 
- .107059E+02 
- .438665E+03 
- .392255E+03 


SHEAR/MOMENT  RATIO  MAGNITUDES  (X,Y,S23MAX/S22MAX,S13MAX/S11MAX) 
.25446E-02  .25446E-02  0.167694E-02  0.113155E-02 

.25446E-02  .97455E-01  0.658002E-01  0.110972E-02 

.97455E-01  .25446E-02  0.166019E-02  0.438142E-01 

.97455E-01  .97455E-01  0.651023E-01  0.428504E-01 

MAX  (MAGNITUDE)  TRANS  NORMAL  STRESS  (X,Y, ESTIMATE) 

.25446E-02  .25446E-02  -.612563E+04 

.25446E-02  .97455E-01  -.579547E+04 

.97455E-01  . 25446E-02  -.588080E+04 

.97455E-01  .97455E-01  -.555064E+04 
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